#####################################################################################
####Name of File: ifb_statespecific.R
####Input dataset: "Unrestricted_SWU.Rdata"
####Output : state specific regression output
####Date: 4/24/2011
####Purpose: regress ifb on jsy and covariates for each state
#####################################################################################



##Run regression on state specific outcome on jsy uptake 
library(Zelig)
library(survey)
library(tcltk)
library(cem)
library(MatchIt)


##set seed
set.seed(02138)



###load data
load("N:/Repdata/Table1/Unrestricted_SWU.RData")


summary(Unrestricted.swu)  #zero NAs in dataset so can go ahead with matching
dim(Unrestricted.swu)  #

Unrestricted.swu$religion2<-NULL

Unrestricted.swu$wealthgroup<-ifelse(Unrestricted.swu$dlhs3_country_quintile<=2,1,
						ifelse(Unrestricted.swu$dlhs3_country_quintile==3, 2,3))

Unrestricted.swu$castegroup<-ifelse(Unrestricted.swu$caste<=2,1,
						ifelse(Unrestricted.swu$caste==3| Unrestricted.swu$caste==4, 2,99))
						
Unrestricted.swu$educgroup<-ifelse(Unrestricted.swu$mateduc_cat<=1,1,
				ifelse(Unrestricted.swu$mateduc_cat==2| Unrestricted.swu$mateduc_cat==3, 2,99))

Unrestricted.swu$agegroup<-ifelse(Unrestricted.swu$matage_num<=2,1,2)

Unrestricted.swu$nbgroup<-ifelse(Unrestricted.swu$nb_cat_num==1,1,
				ifelse(Unrestricted.swu$nb_cat_num==2,2,
				ifelse(Unrestricted.swu$nb_cat_num==3| Unrestricted.swu$nb_cat_num==4,3,99)))

###prepare unrestricted dataset for matching and logit

Unrestricted.swu$mateduc_cat<-ifelse(Unrestricted.swu$mateduc_cat==99,NA, Unrestricted.swu$mateduc_cat)
Unrestricted.swu<-na.omit(Unrestricted.swu)
dim(Unrestricted.swu)  #177603


####State specific regressions - matching


countries.all<-sort((unique(Unrestricted.swu$state)))
all.data<-matrix(list(),nrow=length(countries.all),ncol=1)


for (i in 1:length(countries.all)){

	x<-countries.all[i]
	print(x)
	z<-Unrestricted.swu[Unrestricted.swu$state==x,]
	
	match.output<-matchit(jsy~bpl_card + urban + as.factor(wealthgroup)+as.factor(castegroup)+ as.factor(educgroup)+as.factor(nbgroup)+as.factor(agegroup), data=z, method="exact")
	
	a<-summary(match.output)
	print(a$nn)
	all.data[[i,1]]<-match.data(match.output, weights="weights")

}

#append all datasets together
all.data2<-all.data

###check dimensions and jsy distribution among the states, to look for sample size issues
for(i in 1:length(countries.all)){
print(countries.all[i])
print(dim(all.data2[[i,1]]))
print(table(all.data2[[i,1]]$jsy))
}




#########################################################################################################################################################
###############################IFB



ifb.output<-matrix(list(),nrow=length(countries.all),ncol=1)

firstdiff.ifb<-matrix(nrow=length(countries.all), ncol=4)


##tjese are the states for which the setx runs; the NAs are because there are errors in compatible arguments


listcount<-c(2:5,7,10,12:17,20:25, 28:30, 32:34)


for(i in listcount){

DF1<-as.data.frame(all.data2[[i,1]])

DF1$nb_cat<-as.factor(DF1$nb_cat_num)
DF1$mateduc_cat<-as.factor(DF1$mateduc_cat)
DF1$dlhs3_country_decile<-as.factor(DF1$dlhs3_country_decile)
DF1$caste1<-as.factor(DF1$caste)
DF1$state1<-as.factor(DF1$state)
DF1$dist.char<-as.character(DF1$DIST_mean_hh_p_income3)
DF1$district<-as.factor(DF1$dist.char)
DF1$religion2<-as.factor(DF1$religion.char)

print(unique(DF1$state1))

ifb.logit1<-zelig(ifb~jsy+matage_cat+nb_cat+birth_interval+mult_birth+mateduc_cat+dlhs3_country_decile+caste1+ religion2+res_cat+ 
+district, model="logit", weights="weights", data=DF1, save.data=T)

ifb.output[[1,1]]<-round(exp(coef(ifb.logit1)),2)

x.treat <- setx(ifb.logit1, jsy=1)
x.control <- setx(ifb.logit1, jsy = 0)
s.out <- sim(ifb.logit1, x = x.control, x1 = x.treat, weights=~swu)

firstdiff.ifb[i,1]<-summary(s.out)$qi.stats$fd[1]*100
firstdiff.ifb[i,2]<-summary(s.out)$qi.stats$fd[2]*100
firstdiff.ifb[i,3]<-summary(s.out)$qi.stats$fd[3]*100
firstdiff.ifb[i,4]<-summary(s.out)$qi.stats$fd[4]*100

}



###########punjab needs special run to get rid of 0 level factor in birth-interval
DF1<-as.data.frame(all.data2[[27,1]])

DF1$nb_cat<-as.factor(DF1$nb_cat_num)
DF1$mateduc_cat<-as.factor(DF1$mateduc_cat)
DF1$dlhs3_country_decile<-as.factor(DF1$dlhs3_country_decile)
DF1$caste1<-as.factor(DF1$caste)
DF1$state1<-as.factor(DF1$state)
DF1$dist.char<-as.character(DF1$DIST_mean_hh_p_income3)
DF1$district<-as.factor(DF1$dist.char)
DF1$religion2<-as.factor(DF1$religion.char)

print(unique(DF1$state1))

DF1$birth_interval <- DF1$birth_interval[DF1$birth_interval!="0-11 Months"]
DF1$birth_interval<- factor(DF1$birth_interval)
table(DF1$birth_interval)

ifb.logit1<-zelig(ifb~jsy+matage_cat+nb_cat+birth_interval+mult_birth+mateduc_cat+dlhs3_country_decile+caste1+ religion2+res_cat+ 
+district, model="logit", weights="weights", data=DF1, save.data=T)

ifb.output[[27,1]]<-round(exp(coef(ifb.logit1)),2)

x.treat <- setx(ifb.logit1, jsy=1)
x.control <- setx(ifb.logit1, jsy = 0)
s.out <- sim(ifb.logit1, x = x.control, x1 = x.treat, weights=~swu)

firstdiff.ifb[27,1]<-summary(s.out)$qi.stats$fd[1]*100
firstdiff.ifb[27,2]<-summary(s.out)$qi.stats$fd[2]*100
firstdiff.ifb[27,3]<-summary(s.out)$qi.stats$fd[3]*100
firstdiff.ifb[27,4]<-summary(s.out)$qi.stats$fd[4]*100

###########tripura needs special run to get rid of 0 level factor in birth-interval
DF1<-as.data.frame(all.data2[[31,1]])

DF1$nb_cat<-as.factor(DF1$nb_cat_num)
DF1$mateduc_cat<-as.factor(DF1$mateduc_cat)
DF1$dlhs3_country_decile<-as.factor(DF1$dlhs3_country_decile)
DF1$caste1<-as.factor(DF1$caste)
DF1$state1<-as.factor(DF1$state)
DF1$dist.char<-as.character(DF1$DIST_mean_hh_p_income3)
DF1$district<-as.factor(DF1$dist.char)
DF1$religion2<-as.factor(DF1$religion.char)

print(unique(DF1$state1))

DF1$birth_interval <- DF1$birth_interval[DF1$birth_interval!="0-11 Months"]
DF1$birth_interval<- factor(DF1$birth_interval)
table(DF1$birth_interval)

ifb.logit1<-zelig(ifb~jsy+matage_cat+nb_cat+birth_interval+mult_birth+mateduc_cat+dlhs3_country_decile+caste1+ religion2+res_cat+ 
+district, model="logit", weights="weights", data=DF1, save.data=T)

ifb.output[[31,1]]<-round(exp(coef(ifb.logit1)),2)

x.treat <- setx(ifb.logit1, jsy=1)
x.control <- setx(ifb.logit1, jsy = 0)
s.out <- sim(ifb.logit1, x = x.control, x1 = x.treat, weights=~swu)

firstdiff.ifb[31,1]<-summary(s.out)$qi.stats$fd[1]*100
firstdiff.ifb[31,2]<-summary(s.out)$qi.stats$fd[2]*100
firstdiff.ifb[31,3]<-summary(s.out)$qi.stats$fd[3]*100
firstdiff.ifb[31,4]<-summary(s.out)$qi.stats$fd[4]*100


#########pull all first differences together


firstdiff2.ifb<-as.matrix(round(firstdiff.ifb,4))
rownames(firstdiff2.ifb)<-countries.all

firstdiff.final.ifb<-firstdiff2.ifb[,-2]


#############do the plotting
###create axis names 

##take out NAs for states that couldn't calculate first differences


plotdata<-as.matrix(firstdiff.final.ifb[-c(1,6,8,9,11,18,19,26),])

plotdata2<-plotdata[sort.list(plotdata[,1]), ]

#set up plot area
par(mar=c(7, 8, 3, 2) + 0.1)

##plot age
y1<-rep(1,2)
##start with one state
plot(plotdata2[1,1],1, ylim=c(0,26), xlim=c(-5,78),axes=F, xlab= "In-Facility Birth State Level Treatment Effect (%)")
points(plotdata2[1,2:3],y1, type="o", pch="|")



##draw axes
box()
#axis(1, at=c(2,8), lab=c("IFB"))
axis(2, at=c(1:26), lab=rownames(plotdata2),las=1)
axis(1, seq(-4, 78, by=4), cex=.6)

#Continue for all other states
for(i in 2:26){
y<-rep(i,2)
points(plotdata2[i,1],i, pch="o")
points(plotdata2[i,2:3],y, type="o", pch="|")

}

abline(v=47.04, col="blue")
abline(v=45.51, col="light blue", lty=2)
abline(v=48.41, col="light blue", lty=2)
abline(v=0, col="gray", lty=2)

text(60, 10, "Blue line is national level effect\n with 95% confidence intervals")








#######################################################################################################################################
###########################################################BALANCE CHECKING FOR SELECT STATES#############################################


########################################################################################
###########pre balifbe
Bihar.pre<-Unrestricted.swu[Unrestricted.swu$state=="Bihar",]
head(Bihar.pre)
dim(Bihar.pre)#17517

balifbe1 <- imbalifbe(group= Bihar.pre$jsy, data=Bihar.pre, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.856, LCS=5.6%


###########post balifbe
Bihar<-as.data.frame(all.data2[[5,1]])
head(Bihar)
dim(Bihar)#16822

balifbe1 <- imbalifbe(group= Bihar$jsy, data=Bihar, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.851, LCS=5.9%

########################################################################################
###########pre balifbe
Assam.pre<-Unrestricted.swu[Unrestricted.swu$state=="Assam",]
head(Assam.pre)
dim(Assam.pre)#7612

balifbe1 <- imbalifbe(group= Assam.pre$jsy, data=Assam.pre, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.874, LCS=7.3%

###########post balifbe
Assam<-as.data.frame(all.data2[[4,1]])
head(Assam)
dim(Assam)#16822

balifbe1 <- imbalifbe(group= Assam$jsy, data=Assam, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.874, LCS=7.5%

########################################################################################

###########pre balifbe
Madhya.pre<-Unrestricted.swu[Unrestricted.swu$state=="Madhya Pradesh",]
head(Madhya.pre)
dim(Madhya.pre)#13624

balifbe1 <- imbalifbe(group= Madhya.pre$jsy, data=Madhya.pre, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.845, LCS=8.6%

###########post balifbe
Madhya<-as.data.frame(all.data2[[18,1]])
head(Madhya)
dim(Madhya)#13461

balifbe1 <- imbalifbe(group= Madhya$jsy, data=Madhya, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.844, LCS=8.7%

########################################################################################

###########pre balifbe
Rajasthan.pre<-Unrestricted.swu[Unrestricted.swu$state=="Rajasthan",]
head(Rajasthan.pre)
dim(Rajasthan.pre)#11004

balifbe1 <- imbalifbe(group= Rajasthan.pre$jsy, data=Rajasthan.pre, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.857, LCS=8.1%

###########post balifbe
Rajasthan<-as.data.frame(all.data2[[26,1]])
head(Rajasthan)
dim(Rajasthan)#10793

balifbe1 <- imbalifbe(group= Rajasthan$jsy, data=Rajasthan, 
drop=c("jsy","obs","ifb3","ifb","enm","pnm","ab","nmsb","lnm","nnm","birth_interval","mult_birth","dwu",
"swu","state_cat_num","matage_num","bpl_card_num","nb_cat_num","wealthgroup","castegroup","educgroup","agegroup","nbgroup","weights","subclass","state_cat",
"dweight","dm_coef","sm_coef","cmcage","state","dlhs3_country_quintile","bpl_card"))
balifbe1 
###L1=.855, LCS=8.2%

######################################################################################################################################################


###models that did not run: 2,3,10,11,13,22,24,28,30,38

full.output<-as.data.frame(ifb.output2[1,1])
for (i in 2:length(ifb.output2)){
	full.output<-cbind(full.output,as.data.frame(ifb.output2[i,1]))
}

for (i in 

####haven't tried this
##TABLE 3 estimates

x.treat <- setx(ifb.logit1, jsy=1)
x.control <- setx(ifb.logit1, jsy = 0)
s.out <- sim(ifb.logit1, x = x.control, x1 = x.treat, weights=~swu)

summary(s.out)$qi.stats$fd[1]

firstdiff<-matrix(nrow=












